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ABSTRACT 

Current mesh reduction techniques, while numerous, all pri- 
marily reduce mesh size by successive element deletion (e.g. 
edge collapses) with the goal of geometric and topological 
feature preservation. The choice of geometric error used to 
guide the reduction process is chosen independent of the 
function the end user aims to calculate, analyze, or adap- 
tively refine. In this paper, we argue that such a decou- 
pling of structure from function modeling is often unwise 
as small changes in geometry may cause large changes in 
the associated function. A stable approach to mesh decima- 
tion, therefore, ought to be guided primarily by an analysis 
of functional sensitivity, a property dependent on both the 
particular application and the equations used for computa- 
tion (e.g. integrals, derivatives, or integral/partial differen- 
tial equations). We present a methodology to elucidate the 
geometric sensitivity of functionals via two major functional 
discretization techniques: Galerkin finite element and dis- 
crete exterior calculus. A number of examples are given to 
illustrate the methodology and provide numerical examples 
to further substantiate our choices. 

1. INTRODUCTION 

For function computations carried out on large meshes, mesh 
decimation is an essential first step. Mesh decimation tech- 
niques are distinguished by the cost function they attempt 
to minimize as they collapse edges in the mesh. In this 
paper, we show that given a particular partial differential 
equation (PDE) problem, an analysis of the geometric sen- 
sitivity of the functions involved should guide the choice of 
cost function for pre-computation decimation. 

We consider such function-guided decimation in two realms: 
adaptive finite element methods (AFEM) and discrete exte- 
rior calculus (DEC) methods. There are three main steps in 
such methods: formulating a weak version of the governing 
PDEs, discretizing the problem, and solving the resulting 
linear system. Each step introduces a different type of er- 
ror to the process. Formulating a weak problem may create 
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what is known as model error. Reducing to a linear sys- 
tem causes discretization error. Implementing the numerical 
method inevitably causes some solver error. 

Adaptive finite element methods aim to control solver error 
by selective local refinement of the input mesh (/i-adaptivity), 
the degree of the basis functions (p-adaptivity) , or both (hp- 
adaptivity). While each flavor of adaptive method has met 
success in particular applications, we will show that their ap- 
plicability does not immediately transfer to problems that 
require mesh decimation as a pre-processing step. Mesh dec- 
imation causes a certain loss of geometric information while 
adaptive refinement is an approximation of missing function 
information. For this reason, it is important that the loss 
of function information accrued during mesh decimation be 
bounded a priori so that the adaptive method can have a 
hope of converging to a meaningful result. 

Discrete exterior calculus methods control solver error by 
discretizing the functions and operators of the PDE with re- 
spect to their algebraic relationships. This type of analysis 
leads to specific conclusions about where values of the load 
data and solution data should be assigned or computed; in 
many cases, values belong most naturally somewhere other 
than mesh vertices, e.g. on mesh edges or at the circumcen- 
ters of triangles. Therefore, an error bound on function in- 
formation loss for mesh decimation prior to a DEC method 
must, by necessity, take into account the locations of the 
samples of the various variables in the problem. 

We describe a framework for selecting an appropriate mesh 
decimation technique given a PDE and an approach to solv- 
ing it. The discretization from the AFEM or DEC method 
yields a linear system of the form Ax = b whose solution re- 
quires inverting the matrix A. Our contention is that mesh 
decimation should be guided by an attempt to avoid large 
entries in the matrix A which can make A ill-conditioned 
and hence destabilize the numerical method. 

In Section [2] we discuss prior work on AFEM, DEC, and 
mesh decimation. In Section[3] we first give a general overview 
of AFEM and DEC methods and then explain how each can 
suggest a mesh decimation technique through a variety of 
examples. In Section [4] we describe the cost functions as- 
sociated to two existing techniques as well two novel cost 
functions for use in molecular solvation energetics computa- 
tions. In Section [5] we present initial experimental results 
comparing our technique to prior ones. 



2. PRIOR WORK 

We discuss the three main topics of prior work related to our 
approach: adaptive finite element methods, discrete exterior 
calculus methods, and mesh decimation methods. 

Finite element methods (FEM) have witnessed an explosive 
growth both in the literature and in industrial application 
in the past few decades. Adaptive methods [2] have gained 
traction for their ability to increase local accuracy in a so- 
lution. Beginning with a coarse mesh, AFEM refine by sub- 
dividing certain elements into smaller pieces (/i-adaptivity) 
[T6] , increasing the degree of polynomial approximation on 
some elements (p-adaptivity) [5], or a combination of the 
two (/ip-adaptivity) [12]. Recently AFEM have been ap- 
plied to computational biology disciplines. Baker et al. have 
worked on a parallel implementation of an AFEM to solve 
the Poisson-Boltzmann equation 8 which generates a mesh 
of the molecular surface via a subdivision scheme. Recent 
work by Chen et al. [l2] provides a FEM for the nonlin- 
ear Poisson-Boltzmann equation with rigorous convergence 
estimates. 

Discrete Exterior Calculus (DEC) is an attempt to create 
from scratch a discrete theory of differential geometry and 
topology whose definitions and theorems mimic their smooth 
counterparts. This theoretical foundation allows the canon- 
ical prescription of a discretization scheme for a given PDE 
problem that enforces topological constraints combinatori- 
ally instead of numerically, thereby providing for increased 
robustness in implementation. This approach has been em- 
ployed by an increasing number of authors in recent years 
to develop multigrid solvers [9], s olve Darcy flow problems 
[22] , and geometrize elasticity [29] . For a complete introduc- 
tion to DEC theory, see Hirani~[2l] and Desbrun et al. [15] . 
In this paper, we give a brief introduction to the theory in 
Section [3] and an example in Section pO] 

Mesh decimation techniques are abundant in geometry pro- 
cessing literature. A useful survey of many methods was 
given by Heckbert and Garland [20] and a more recent book 
by Luebke et al. [27] provides a thorough overview of mesh 
simplification techniques. Bajaj and Schikore 6 have given 
an error bounded mesh decimation technique for 2D scalar 
field data. We focus on mesh decimation of unstructured 
surface meshes via edge contraction including the approaches 
of Garland and Heckbert [l8] and Lindstrom and Turk |23[ 
1 24] . These are explained in Section [4] 

3. METHODOLOGY FOR ELUCIDATING 
GEOMETRY SENSITIVE FUNCTIONALS 

We begin with a general problem: find u £ V such that 

Lu = f on Q, (1) 

where L is a linear operator, V is the appropriate solution 
space for the problem, and Q is a simplicial complex embed- 
ded in R 3 . There are two main techniques used to discretize 
this into a linear system: Galerkin methods and Discrete 
Exterior Calculus (DEC) methods. We describe how each 
could be used and how a mesh decimation technique should 
be chosen accordingly. 

The Galerkin finite element method begins by putting (fTl) 



into the weak form: find u £ V such that 

a(u,v) = (f,v) L 2 VveV, (2) 

where a is the operator L phrased as a bilinear form (usually 
symmetric) and / is treated as a functional (/, -) L 2 on V. 
An appropriate finite dimensional subspace Vh C V is chosen 
and an answer to the following discretized problem is sought: 
find u G Vh such that 

a(u,v) = (f,v) L 2 Vv e V h . 

Since Vh is finite-dimensional, we can fix a basis : 1 < 
i < n} of Vh- The size of the basis is proportional to the 
number of elements in the mesh Q. Write u — X^=i Uj4>j, 
Kij = a{<t>j^i) and F t = (/,&). Set U = (Uj), K = (K tj ) 
F = (Fi). Then solving (pi) over Vh is the same as solving 
the matrix equation 

KTJ = F. (3) 
For a proof and detailed discussion, see [1 1] . 

A significant amount of care goes into the selection of Vh 
to ensure that the method is both well-posed and stable. 
"Well-posed" means the system has a unique solution and 
"stable" means there exists a constant C > independent of 
h such that 

\\u — Uh\\v<C inf \\u — Wh\\v- 

In other words, a method is stable if the error between the 
true solution u and approximate solution Uh is bounded 
above uniformly by a constant multiple of the minimal ap- 
proximation error for Vh- The famous Babuska inf-sup con- 
dition 1 is often used to simultaneously prove both well- 
posedness and stability of a FEM and hence provide an a 
priori bound on solver error. We note, however, that the 
stability error bound does not account for error due to naive 
mesh decimation. 

If the mesh Q has too many elements, (|3| will be too large 
for the solver, making decimation necessary. For decimation 
to be useful, however, it must not create very large or small 
entries in IK which might make IK ill-conditioned. Since the 
entries of IK are a functional a(-, •) on the basis functions 
mesh decimation must be guided by the geometry-sensitive 
components of a as opposed to the geometry of Q alone. 
Such components are necessarily problem-specific as we elu- 
cidate in examples presented in Sections |3.1| and |3.2| 

The method of Discrete Exterior Calculus is an alternative 
approach which focuses on correctly discretizing the opera- 
tor L instead of the solution space V. The viewpoint pro- 
vided by differential geometry and topology reveals how this 
ought to be done. Common operators such as grad, curl, and 
div are all manifestations of the exterior derivative operator 
d in dimensions 1,2, and 3, respectively. Equations relating 
quantities of complementary dimensions, such as the con- 
stitutive relations in Maxwell's equations, involve a Hodge 
star operator * which provides the canonical mapping. The 
Laplacian operator A can be written as Sd + dS where 5 is 
the coderivative operator, defined by S := Each opera- 
tor has a discrete version designed to mimic the properties 
of its smooth counterpart. The discrete versions of the op- 
erators are written as matrices whose entries depend only 
on the topology and geometry of the mesh Q. 



To solve |T]), an analysis is made as to the dimension of 
u as a /c-form based on either the problem context or the 
type of operator acting on it. The div operator in 3D, for 
example, acts on 2-forms while grad acts on 0-forms. The 
variable u is replaced by a vector u with one entry for each 
/c-simplex in the mesh of Q and the operator L is replaced 
by its discrete counterpart, written as a matrix L. The load 
data / is converted to a vector / accordingly. This yields 
the equation 



which can then be solved by linear methods. 



(4) 



Again, it may be necessary to decimate Q so that the lin- 
ear system Q is tractable on a computer. The size of the 
entries of L depend heavily on the geometry-sensitive op- 
erators such as * and S and less on the topology sensitive 
operators such as d. Therefore, to prevent an ill-conditioned 
L, mesh decimation must be guided based on the definition 
of the discrete operators as opposed to the definition of the 
solution space. We discuss an example in Section [33] 



3.1 Poisson-Boltzmann Electrostatics 

The Poisson-Boltzmann equation (PBE) describes the at- 
traction between solvated molecules. We describe its lin- 
earized version according to t he formulation given by Lu, 
Zhang, and McCammon in [26], which is believed to be a suf- 
ficient approximation for electrostatics computations. Let 
Q C R 3 be a domain indicating interior molecular regions 
with point charges qi,. . . ,qN located at n,. . . ,rjv G O. The 
linear PBE is 



V¥ nt (r p ) 



1 

^2 QkS(r p - r fc ), pGO, 



Cint 



V 2 ext (r p ) = ^ 2 ext (r p ), O, 



(5) 



(6) 



where mt and ext are the electrostatic potentials on the 
interior and exterior of Q, £mt is the interior dielectric con- 
stant, 5 is the Dirac distribution, and k is the inverse of the 
Debye-Hiickel screening length. Values for the constants Cmt 
and k are determined experimentally. At the boundary <9Q, 
the surface potential / should satisfy / = mt = ext with 
normal derivative h = <9(/> ext /dn. The equations are put into 
integral form and discretized, reducing the problem to a set 
of linear equations of the form 







(7) 



where Q is the initial data of point charges and locations. 
The four (sub) matrices A, B, C and D have entries 



XdA 



where Et is a facet of Q, t indexes over a small neighborhood 
of mesh elements, and the integrand X is one of the following: 



X G { G(xi,Xj), 



dG(xj,Xj) 
dn 



l\Xi , Xj ) 



du(xi, Xj 
dn 



The functions G and u are the Green functions for |5| and 
([6]), respectively. They are given by 

nt \ 1 A ( \ GXp(-Krij) 

G(xi,Xj) = - and u{Xi,Xj) = 



where nj = \x% — Xj\. Hence, the terms of the submatrices 
in |7| decay like 1/rij at worst. To avoid a blowup of these 
terms, we use a cost function j vh that attempts to preserve 
mesh element quality as a way of avoiding small nj values. 
We describe f p b in Section [43] 



3.2 Generalized Born Electrostatics 

A recent approach by Bajaj and Zhao \7\ computes molec- 
ular solvation energetics and forces by using a Generalized 
Born (GB) model instead of a Poisson Boltzmann model 
of electrostatic solvation. While the PB model begins with 
the PBE, a description of the electrostatic potential over 
the whole domain, the GB model begins with a model of 
the solvation energy of a single atom in a given medium. 
The electrostatic solvation energy of a molecule is defined 
in terms of the pairwise interaction between these atomic 
energies: 



G 



pol 



QiQi 



Here, r — where e v and e w are the solute and solvent 

dielectric constants, qi and Ri are the charge and effective 
Born radius of atom i, and nj is the distance between atom 
i and atom j. The success of a GB method hinges upon effi- 
cient and accurate approximation of the effective Born radii 
Ri. Bajaj and Zhao use the surface integration technique 
given in [19] . This gives the expression 

4tt J v |r-x^| 4 

where Y is the solvent-molecular interface, x» is the center 
of atom z, and n(r) is the unit normal of the surface at 
r. The position vector r ranges over Y. The integral is 
approximated by 



R 



1 N 
4tt ^ 



(rjfc - Xi) • n(rjfe) 



4 dS, i = l,...,M, 



where the are the Gaussian quadrature nodes with weights 
Wk lying on a triangular mesh approximating the surface 
T. Therefore, the computation of R^ 1 will be sensitive to 
changes in the position of Gauss points relative to the near- 
est atomic centers, i.e. changes in the computed values of 
|r — Xi|. Accordingly, we design a cost function f g b to mini- 
mize the cumulative change in |r — x^| values. This function 
is described in Section [4. 4 1 and compared experimentally to 
other cost functions in Section [5] 

3.3 Darcy Flow 

Recent work by Hirani et al. [22] uses a DEC method to 
model Darcy flow, a description of the flow of a viscous fluid 
in a permeable medium. The governing equations under the 
assumption of no external body force are given by 



div/ 
/ 







in 

cj) in Q, 
on dQ, 



where / is the volumetric flux, k > is the coefficient of 
permeability, \i > is the coefficient of viscosity, cj) : Q — > R 
is the prescribed divergence of velocity, and ip : — > R is 
the prescribed normal component of the velocity across the 




Figure 1: Notation for collapse of edge (^1,^2) to 
v. Black vertices and edges are unchanged in the 
collapse while red vertices and dashed edges may 
change. 



boundary. They discretize the operators V and div based 
on DEC theory to arrive at the linear system 
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Here, is the discrete exterior derivative operator that acts 
on /c-cochains and is a diagonal matrix representing the 
Hodge Star operator on /c-cochains. 

We now consider the effect of a priori mesh decimation for 
this scheme. The DEC analysis used to derive this method 
requires that the solution [/ p] T provide values of the flux 
/ on (n — l)-simplicies (i.e. edges in triangle meshes and 
triangles in tetrahedral meshes) and values of the pressure 
p at the circumcenters of n-simplicies. 

For the pressure values to have any meaning, the mesh must 
be well- centered meaning the circumcenter of each simplex 
must lie in the interior of the simplex. Since this criterion 
is often violated by meshing schemes (e.g. an obtuse tri- 
angle is not well-centered), pressure is assigned instead to 
the bary centers of n-simplicies. The authors point out that 
the error introduced by this modification prevents the ex- 
act representation of linear variation of pressure over the 
domain. Therefore, an appropriate cost function for this 
method should be weighted to favor the creation of simpli- 
fied meshes with good quality elements (e.g. elements with 
good aspect ratios). This would minimize the distance be- 
tween the bary center and circumcenter, making the calcula- 
tions more robust. 

4. DESCRIPTION OF COST FUNCTIONS 
4.1 Quadratic Error Cost Function 

The quadric error measure we use comes from Garland and 
Heckbert [18] . First, an error metric A(v) is established 
for each vertex v, based on the planes P(v) that contain 
the triangles incident to v. A plane p £ -P(v) given by 
ax + by + cz + d = is represented as [a b c d] T . The 

v z 11 T . Then the error 



vertex v is represented as [v x 
metric A(v) is defined by 



A(„)=v T £ PP 1 

\pGP(v) 



For a vertex v;, let Q; = ^ p€ p( Vi ) PP 
collapsing edge (vi, v 2 ) to some point v is 



Then the cost of 



/ge(vi,V 2 ;v) 



+ Q 2 )v. 



We use the publicly available software QSlim [17] to im- 
plement this cost function. For a given edge, the program 
attempts to find an optimal placement of v by solving a 
certain linear system derived from the Qi. If the matrix as- 
sociated to this system is not invert ible, it tries to place v 
optimally on (vi, V2). If this fails, it sets v to be either vi, 
V2, or the midpoint of the edge, whichever minimizes f qe . 

4.2 Volumetric Error Cost Function 

The volumetric error measure we use comes from Lindstrom 
and Turk [23[ |24| . First, we establish the notation for the 
signed volume V of a tetrahedron bounded by a vertex v and 
the vertices Vq% v^ of a triangle U. As before, vertices 
are written as four component vectors, e.g. v is represented 
by [v x v y v z 1] T . Then V is defined by 



VXvjVoSvjSva*) = ^det(v v£ v** v£) =: ^Gt.v 

where Gt i is a 1 by 4 matrix defined by the above equation. 
The cost of collapsing edge (vi, V2) to some point v is 



fvoiiy±, v 2 ;v) := -v 



where i indexes over triangles U incident to at least one of 
{vi,V2}. This cost function is ultimately quite similar to 
f qe except that f qe weights the distance between v and a 
plane by the area of the triangle defining the plane while 
f V oi weights it by the square of the triangle area. As is 
explained in |23[ |24| , this weighting better serves the goal 
of volume preservation. We use a package from the publicly 
available software TeraScale Browser [25] to implement this 
cost function. 

4.3 Poisson Boltzmann Cost Function 

In Section 3.1 we discuss how Poisson Boltzmann (PB) com- 



putations are sensitive to Gaussian quadrature points com- 
ing into close proximity. Hence, we define a cost function 
f p b which penalizes such occurrences. 

We fix the following notation for the collapse of edge (vi , V2) 
to the point v. Consider the union of triangles incident to vi 
or v 2 . These are the only vertices, edges, and triangles whose 
geometry may be changed by the edge collapse. Although 
all these objects lie in R 3 , their generic connectivity infor- 
mation is captured in R 2 by Figure [l] (a) . Taking (vi,V2) 
to be vertical with vi at the bottom, the triangle to the left 
(resp. right) of the edge has vl (vj?) as its third vertex. 
We proceed from vz, to vr along the upper (resp. lower) 
vertices labeling them v 2 , v|, . . ., v^ (vi, vf, . . ., vf (D for 
down)). Set v§ := v? := w L and v^ +1 := vf +1 := w R . We 
denote the centers of the triangles not adjacent to (vi,V2) 
as 

1 

3 



c 2 



( V 2 + V 2 + V 2 



U ■ 



0,1, 



,u, 



d 
Cl 



1 d 1 d-\- 
1 + V 1 + Vi 



') 



0,1, 



,D. 



The collapse operation moves vi and V2 to v. The result is 
shown in Figure 0(b). The {vf} and {v£} are unchanged, 
but the new centers are {c%} and {cf} where v replaces V2 
or vi in the expression of the center. Re-index cf and C2 as 



If the triangles indexed by the Ci are of good quality, meaning 
nearer-to-equilateral, their Gaussian quadrature points will 
be better spaced. We approximate the quality of a triangle 
Tby 



l(T) maxa(T) 
s(T) mina(T) ' 



where l(T) (resp. s(T)) denotes the longest (shortest) side 
and maxa(T) (resp. mina(T)) denotes its maximum (minu- 
mum) angle. The best quality triangles have the minimum 
q value of 2. We denote triangles in Figure [I] by their center 
Ci or Ci. The cost function is then defined to be 

/p»(vi,v 2 ;v) := ^>(TcJ - q(T Ci ). 

i 

To further improve element quality, we decimate in stages 
and run a quality improvement code based on geometric flow 
[28] in between stages. 

4.4 Generalized Born Cost Function 



In Section 3.2 we discuss how the Generalized Born (GB) 
computations are sensitive to changes in the location of 
Gaussian quadrature points {d} relative to the atomic cen- 
ters {xj} of the molecule in question. Hence, we define a cost 
function f g b which penalizes edge collapses with a larger cu- 
mulative change in \a — Xj\ values. Since it would be too 
computationally expensive to compute the complete change 
for every edge collapse, we use a restricted set of pertinent 
{ci} and {xj} values described below. 



4.3 



as our set of per- 



We take the {d} described in Section < 
tinent Gaussian quadrature points. We set {xj} to be those 
atomic centers lying within a fixed distance p of either v\ or 
V2. Since Born radii are on the order of 1-2 A, p should be 
set between 2 and 5 to capture a manageable, non-empty set 
of nearby atoms. In the future, we will devise a parameter 
sweep to optimize the value of p. 

We want to minimize the atomic center functional f ac 



£1 



Ci Xj 



which can be re-written as 



E 



T 

Ci Ci 



■ cj Ci + 2(cJ 



Ci )Xj 



All the variables in the above expression are known, save 
for the Ci which are linear functions of v. We define the 
GB-dependent cost function to be 

f 9 b(vi, v 2 ; v) := |vi - v 2 | + A/ ac (vi, v 2 ; v). 

The first term is used to promote the collapse of shorter 
edges and thereby improve triangle quality. The weight fac- 
tor A is chosen so that the two terms are of the same order 
of magnitude. 




Figure 2: Top: Surface rendering of mAChE before 
decimation. The original mesh of the surface has 
about 200,000 vertices and 400,000 faces. In the 
inset, the fine mesh is visible. This pocket region of 
the molecule aids in its biological function. Bottom: 
The mesh of the pocket after decimation to 25,500 
faces using f qe (left) and a modified version of f g b 
(right). 



5. EXPERIMENTAL RESULTS AND CON- 
CLUSIONS 

To test the validity of our claims, we work with Mouse 
Acetylcholinesterase (mAChE). This macromolecule serves 
an important regulatory function as it terminates the ac- 
tion of the neurotransmitter acetylcholine (ACh). The ini- 
tial mesh of the molecular surface is generated from a Pro- 
tein Data Bank (PDB) [10] fi le of the molecule using our 
in-house software TexMol |13|. We show a picture of the 
initial mesh in Figure [2] We decimate this mesh using the 
cost function f qe and a modified version of f g b which uses f qe 
instead of Jy 1 — V2 1 . We use a Dynamic Packing Grid data 
structure p\ to efficiently compute the set {xj} of nearby 
centers for each mesh edge. We set p = 5 and A = 10 ~ 8 and 
compute the polarized and non-polarized energy for each 
mesh using the nFFGB code described in 7 . The meshes 
are only marginally different as shown in Figure [2] and thus 
produce similar energy values as shown in the charts in Fig- 
ure [3] With further experimentation and parameter sweeps, 
we believe f g b will begin to out-perform f qe . Still, Figure [3] 
shows that f vo i is a decidedly worse choice for non-polarized 
energy computations as it does not respect the functional 
sensitivity of the problem. In future work, we will also im- 
plement f p b and use PB-CFM code by Bajaj and Chen [2] 
to compute and compare PB energetics. 
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Figure 3: Effect of cost function used for decimation 
on computed GB energy values for mAChE. 
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